gampmatlabfandomcom-20200213-history
Example: Estimation of a Gaussian Vector
Main -> Users' Guide and Examples -> Estimation of a Gaussian Vector. Next: Estimating a sparse vector. Program The code for the example is in examples/basic/simpleAWGN.m in the directory gampmatlab/trunk/code. To run the code, simply cd to that directory and run simpleAWGN at the MATLAB prompty. You will see a scatter plot of the components of a vector against two estimates: the optimal least-squares (LS) estimate and the GAMP estimate. In addition, the mean squared error (MSE) for both methods will be printed. Problem Description In this first problem, we will show how to use the gampEst function to estimate a Gaussian random vector \mathbf{x} from random linear measurements of the form \quad \mathbf{y}=\mathbf{z} + \mathbf{w}, \quad \mathbf{z}= \mathbf{A}\mathbf{x}, where \mathbf{A} is a Gaussian random matrix and \mathbf{w} is Gaussian noise. The solution to this problem is, of course, standard and be performed easily without GAMP: For example, if the components of \mathbf{x} and \mathbf{w} are i.i.d., Gaussian zero-mean with variances \sigma^2_x and \sigma^2_w , respectively, then the MMSE estimate is given by the least squares solution: \quad \widehat{\mathbf{x}}_{\rm LS} = \sigma^2_x\left( \sigma^2_x\mathbf{A}'\mathbf{A} + \sigma^2_wI\right)^{-1}\mathbf{A}'\mathbf{y}. \quad (1) Nevertheless, we will show how to get an approximate solution to this problem using GAMP simply to illustrate the syntax and usage of the MATLAB command. MATLAB Program Description Generation of the Data At the top of the program simpleAWGN.m, you will see the lines: % Set path to the main directory addpath('../../main/'); % Parameters nx = 1000; % Number of input components (dimension of x) nz = 2000; % Number of output components (dimension of y) snr = 20; % SNR in dB. The first line adds the necessary directories to the MATLAB path. You can also set the MATLAB path permanently if you intend to use GAMP significantly. The remaining lines set the parameters of the problem. Feel free to change these to experiment. The next lines generate a random Gaussian vector \mathbf{x} : % Create a random Gaussian vector xmean0 = 0; xvar0 = 1; x = normrnd(xmean0, sqrt(xvar0),nx,1); The matrix \mathbf{A} and the noise vector \mathbf{w} are generated with Gaussian i.i.d. components, with the noise scaled according to the parameter snr: % Create a random measurement matrix A = (1/sqrt(nx))*randn(nz,nx); % Compute the noise level based on the specified SNR. Since the components % of A have power 1/nx, the E|y(i)|^2 = E|x(j)|^2 = xvar. wvar = 10^(-0.1*snr)*xvar0; % Generate the noise w = normrnd(0, sqrt(wvar), nz, 1); y = A*x + w; Solution via GAMP The main GAMP command to perform the estimation of \mathbf{x} from the measurement vector \mathbf{y} appears near the bottom of the file: xvar = gampEst(inputEst, outputEst, A, opt); The gampEst function takes as input arguments: * inputEst is a class of type EstimIn that describes the prior p(x_j) on the components of the unknown vector \mathbf{x} ; * A is a matrix for the transform for the relation \mathbf{z}=\mathbf{A}\mathbf{x} . This can also be specified by a LinTrans class as we will see in subsequent examples. * outputEst is class of type EstimOut describing the transition probabilities for the measurements p(y_i|z_i) ; and * opt is a set of options. The outputs xhat and xvar are the estimates of the mean and variance of the vector \mathbf{x} . You can, in fact, get other outputs to estimate the full marginal posteriors of each the variables -- we will illustrate this in a later example. The idea of specifying the prior and measurement channel via the MATLAB classes inputEst and outputEst may be confusing at first. However, as we will see in this sequence of examples, the use of classes allows the algorithm to be very general, modular and capable of handling many different types of distributions. If you are not familiar with MATLAB classes or object-oriented programming, you can read more on the MATLAB object-oriented webpage. However, the gampmatlab package contains pre-built classes for many common input and output distributions. For most problems, you can just use one of these pre-built classes and not have to build one from scratch. In this particular example, since the input is a Gaussian random vector, the class inputEst is created using the pre-built AwgnEstimIn class with the line: inputEst = AwgnEstimIn(xmean0, xvar0); The output estimator class is also created via a pre-built class. Specifically, the line outputEst = AwgnEstimOut(y, wvar); creates outputEst class using the function AwgnEstimOut for an AWGN measurement channel. Note that the observation vector y is passed to the estimator class -- not the gampEst function. With the estimator classes defined, the program then calls gampEst which returns an estimate vector xhat. The GAMP estimate is compared to the MMSE estimator. For this problem, the MMSE estimate is given by the solution in equation (1) above to the least-squares problem. This estimate is computed in the MATLAB code with the lines: xhatLS = xvar0*(A'*A*xvar0 + wvar*eye(nx))\(A'*y); The program concludes by plotting x, the GAMP estimate xhat and LS estimate xhatLS. The MSE for the GAMP and LS estimates are also printed. You should see that GAMP estimate obtains an error close to the LS estimate.